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Abstract 

Background: Osteoporosis affects 200 million people worldwide and places an enormous economic burden on 
society. We aim to identify the feature genes that are related to osteoprotegerin in osteoporosis and to perform 
function analysis with DNA microarray from human bone marrow. 

Methods: We downloaded the gene expression profile GSE35957 from Gene Expression Omnibus database 
including nine gene chips from bone marrow mesenchymal stem cells of five osteoporotic and four non- 
osteoporotic subjects. The differentially expressed genes between normal and disease samples were identified by 
LIMMA package in R language. The interactions among the osteoprotegerin gene (OPG) and differentially expressed 
genes were searched and visualized by Cytoscape. MCODE and Bingo were used to perform module analysis. 
Finally, GENECODIS was used to obtain enriched pathways of genes in an interaction network. 

Results: A total of 656 genes were identified as differentially expressed genes between osteoporotic and 
non-osteoporotic samples. IL17RC, COL1A1, and ESR1 were identified to interact with OPG directly from the 
protein-protein interaction network. A module containing ERS1 was screened out, and this module was most 
significantly enriched in organ development. Pathway enrichment analysis suggested genes in the interaction 
network were related to focal adhesion. 

Conclusions: The expression pattern of IL17RC, COL1A1, and ESR1 can be useful in osteoporosis detection, which 
may help in identifying those populations at high risk for osteoporosis, and in directing treatment of osteoporosis. 
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Background 

Osteoporosis is a chronic disease involving multiple fac- 
tors, and the incidence of osteoporosis in senile people 
and postmenopausal women is rising as the population 
ages, thereby adding to societal problems of health [1]. 
Osteoporosis affects 200 million people worldwide. 
Among those affected, approximately 80% are women 
aged 60 years or older [2]. Current treatment of osteo- 
porosis is mainly by drugs, but it is high cost, time con- 
suming, the drugs have many side effects, and the 
curative effect is not ideal In recent years, more and 
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more scientists have committed to stimulate stem cells 
to differentiate into osteoblasts for the treatment of 
osteoporosis. 

Mesenchymal stem cells have the potential to differenti- 
ate into the osteoblasts [3] or repair bone tissue by 
the lineage or chondrocyte differentiation method [4]. 
Osteoprotegerin (OPG), a member of the tumor necrosis 
factor (TNF) receptor family, suppresses the coupled 
process of skeletal turnover. OPG functions as a decoy re- 
ceptor for osteoclast differentiation factor or as a receptor 
activator of nuclear factor kB (RANK) ligand [5]. Osteo- 
clast differentiation factor promotes bone resorption by 
enhancing the formation and activation of osteoclasts 
when it binds to RANK on hematopoietic osteoclast pro- 
genitor cells as well as on mature osteoclasts [6]. 
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At present, high-throughput screening of differentially 
expressed genes and function identification can get the ex- 
pression profile of bone marrow mesenchymal stem cells, 
and find its differentiation mechanism. But the study of 
high-throughput screening is rare due to expensive equip- 
ment and annotation probe. In this paper, based on a 
group of bone marrow mesenchymal stem cells in gene 
expression profile data, we studied marker genes closely 
linked with osteoprotegerin of osteoporosis in hopes of 
being able to treat osteoporosis through osteoprotegerin 
in bone marrow mesenchymal stem cells. 

Methods 

Affymetrix microarray 

GSE35957 was downloaded from Gene Expression 
Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/ 
geo/), which is based on GPL570 [HG-U133_Plus_2] 
Affymetrix Human information Genome U133 Plus 2.0 
Array Platform (Affymetrix, Santa Clara, CA, USA). 
Microarray probe annotation information was down- 
loaded from the Affymetrix Company, including all 
AffymetrixATHl(25K) gene chip probe information, and 
the probe annotation information files of the platform. 
A total of nine gene chips from mesenchymal cell sam- 
ples, including five gene chips from osteoporosis pa- 
tients and four gene chips from non-osteoporosis 
samples, were used for analysis. 

Data preprocessing and analysis of differentially 
expressed genes 

The original data were preprocessed by Affymetrix [7,8] 
package in R language. LIMMA [9] package in R language 
was used to identify the differentially expressed genes be- 
tween the expression profile of five osteoporosis patients 
and four non-osteoporosis samples. Multiple testing cor- 
rection was performed by Bayesian method [10]. An FDR 
<0.01 and | logFC | >1 were chosen as thresholds for 
screening the differentially expressed genes. 

Prediction of interaction between differentially expressed 
genes 

Differentially expressed genes play a role through inter- 
acting with each other. Therefore, we used HitPredict 
software (http://hintdb.hgc.jp/htp/) to search the differ- 
entially expressed genes that can interact with OPG 
gene. HitPredict is a resource for high confidence 
protein-protein interactions. It collects protein-protein 
interactions from IntAct, BIOGRID and HPRD data- 
bases; annotates these interactions; and assigns a reliabil- 
ity score for each interaction according to the likelihood 
ratio using naive Bayesian networks combining se- 
quence, structure and function annotations of the 
interacting proteins [11]. So far, HitPredict has 239584 
protein-protein interactions across nine species, 168458 



of which are predicted to be of high confidence. This 
study used the protein-protein interactions with high 
confidence to find interactions between the differentially 
expressed genes, and used the Cytoscape [12] to 
visualize the interaction relationships. 

Module analysis of interaction network 

MCODE (Molecular Complex Detection) detects 
densely connected regions in large protein-protein inter- 
action networks that may represent molecular com- 
plexes. In this study, we used MCODE to mine the 
modules from the protein-protein interaction network 
with degree >2. Further, we used Bingo [13] to annotate 
each module based on the hypergeometric distribution 
(FDR <0.05). 

Pathway enrichment analysis of interaction network 

GENECODIS was used to perform biological pathway 
enrichment analysis of all genes in the interaction network 
with FDR <0.05. GENECODIS is a function analysis tool 
of gene, and it integrates different information resources 
(GO, KEGG or SwissProt), searches and arranges gene set 
annotation by statistical significance [14]. 

Results 

Screening differentially expressed genes 

The original data were preprocessed by Affymetrix pack- 
age in R language to remove systematic bias. The results 
show a black line in each box at approximately the same 
level, which indicates an excellent degree of standardi- 
zation (Figure 1). After preprocessing, the normalized 
expression profile data were differentially compared, and 
656 differentially expressed genes exceeding the differ- 
ence threshold (FDR<0.01and |logFC|>l) were screened 
out, including 71 downregulated genes and 585 
upregulated genes. 

Prediction of interaction between OPG and differentially 
expressed genes 

The software HitPredict was used to search all the differ- 
entially expressed genes that interacted with the OPG 
gene. There were 485 interactions among the 656 differ- 
entially expressed genes. Three differentially expressed 
genes (IL17RC, COL1A1, and ESR1) were found to inter- 
act with the OPG gene directly (TNFRSF11B) (Figure 2). 

Module analysis of interaction network 

MCODE was used to mine the densely connected modules 
in protein-protein interaction network. With degree >2, 
one module was screened out (Figure 3). This module con- 
tains 16 nodes, and one of them is ESR1. By using BINGO 
to perform functional annotation, 16 significant GO terms 
were screened out, of which G048513 (involved in organ 
development) is the most significant (Table 1). IGF1R, 
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Figure 1 Expression data after the standardization. The gray box represents the four normal human bone marrow mesenchymal stem cell 
samples, and the white box represents five osteoporosis samples. The black line in each box is the median of data, and its distribution can 
determine the standardization degree of the data. When the black lines are on approximately the same level, this indicates an excellent degree 
of standardization. 




Figure 2 Interaction network among differentially expressed genes and the osteoprotegerin gene (OPG). The gray circular node 
represents the differentially expressed gene, the black triangular node represents the OPG gene and the black square node represents the gene 
directly interacting with OPG. 
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Figure 3 Module diagram of ESR1 gene. The square box represents the ESR1 gene (directly interacting with the osteoprotegerin gene {OPG)), 
and the other circular nodes represent differentially expressed genes. 



APP, BGN, BMP1, ESR1, LOX, ADAMTS2, and TGFB1 
genes were enriched in G048513. 

Pathway enrichment analysis of interaction network 

GENECODIS was used to perform pathway enrichment 
analysis of all differentially expressed genes in the inter- 
action network, and five pathways were significantly 
enriched (Table 2). The pathway of hsa04510: focal ad- 
hesion, involving 21 differentially expressed genes was 
the most significant (FDR = 2.78E-09). The other signifi- 
cantly enriched pathways were the neurotrophin signal- 
ing pathway, regulation of actin cytoskeleton, pathways 
in cancer and the MAPK signaling pathway. 

Discussion 

By the comparison of gene chips from five osteoporosis 
patients and four normal samples of bone marrow stem 
cell, we identified genes (IL17RC, COL1A1, and ESR1) 
that directly interact with the OPG gene. 

The IL17RC (the interleukin 17 receptor C) gene is a 
growth factor that encodes a single-pass type I mem- 
brane protein as extracellular antagonists to cytokine 
signaling [15]. IL-17s and their receptors (IL17RC) pro- 
duced in response to compressive force may affect 
osteoclastogenesis through the expression of RANKL 
and OPG [16]. IL-17RC can also promote bone and joint 
damage through induction of matrix metalloproteinases 
and osteoclasts, and stimulate osteoclastic resorption 
through osteoblasts by inducing receptor activator of 



nuclear factor kB ligand (RANKL) expression [17]. The 
mechanism by which IL17RC interacts with OPG needs 
further study. 

COL1A1 (collagen Type I) is a constituent of the extra 
cellular matrix in connective tissue of bone, skin, ten- 
don, ligament and dentine. It is mostly produced and se- 
creted by osteoblasts and fibroblasts. Mutations in this 
gene are associated with osteogenesis imperfecta types I 
to IV, idiopathic osteoporosis and Caffey Disease [18]. In 
a population-based sample of 1,778 postmenopausal 
women, COLIA1 genotypes of G/G homozygotes (SS), 
G/T homozygotes (Ss), and T/T homozygotes (ss) is as- 
sociated with reduced bone density and predisposes 
women to osteoporotic fractures [19]. COLIA1 Spl(G- 
>T) polymorphism appears to be an important marker 
for low bone mass and vertebral fracture, raising the 
possibility that genotyping at this site may be of value in 
identifying women who are at risk of osteoporosis [20] . 

ESR1 (estrogen receptor 1) localizes to the nucleus 
and plays a role in tissues such as bone, and is involved 
in pathological processes including osteoporosis, endo- 
metrial cancer, and breast cancer. The genes (ESR1, 
BMP1, and IRS1), which are differentially expressed in 
the tibiae of wild type (WT) mice have recognized roles 
in bone metabolism or have been linked previously to 
osteogenesis [21]. Functional annotation showed differ- 
ential expression of ESR1, ESR2, PGR and BGN genes re- 
lated to estrogen metabolism and organ development 
and that these four genes can interact with each other; 
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Table 1 List of module function 



GO-ID 


FDR 


x a 


Description 


Genes in test set 


48513 


1.24E-02 


8 


Organ development 


IGF1R, APP, BGN, BMP!, ESR1, LOKADAMTS2, TGFB1 


10467 


1.24E-02 


7 


Gene expression 


APP, HNRNPUL1, SRSF1 1, PFBP), ESR1,ADAMFS2, HNRNPA 


9888 


1.24E-02 


6 


Tissue development 


IGF1R, APP, BMP], ESR1, ADAMFS2, FGFB1 


16070 


1.34E-02 


6 


RNA metabolic process 


APP, HNRNPUL 1, SRSF11, PFBP1, ESR1, HNRNPA 


9987 


z.Uyb-Uz 


1 6 


Cellular process 


QhADI DTQD1 CDCC7 7 PCD1 UMDMDA TmQI ir~C1D ADD AD1A1 

DlVirl, rlDrl, jnjr / /, Lbnl, HIMnNrA, /(jro/, lur/n, Arr, dO/V, ArzAI, 
HNRNPUL1 ,HGS, LOX, ADAMF52, EPN1, AP2M1 


6807 


2.09E-02 


8 


Nitrogen compound metabolic process 


APP, BGN, HNRNPUL 1, SRSF11, PFBP1, ESR1JGFB1, HNRNPA 


44260 


2.24E-02 


10 


Cellular macromolecule metabolic process 


IGF1R, APP, BGN, HNRNPUL! , SRSF11, PFBP1, E5R1, LOX, FGFB1, HNRNPA 


9653 


2.24E-02 


6 


Anatomical structure morphogenesis 


IGF1R, APP, BGN, BMP!, ESR1, FGFB1 


44238 


2.30E-02 


12 


Primary metabolic process 


IGF1R, APP, BGN, BMP!, HNRNPUL 1, SRSF11, PFBP), ESR1, LOX, ADAMF52, 
FGFB1, HNRNPA0 


6139 


2.30E-02 


7 


Nucleobase, nucleoside, nucleotide and 
nucleic acid metabolic process 


APP, HNRNPUL1, SRSF11, PFBP1, ESR1, FGFB1, HNRNPA0 


48731 


2.38E-02 


8 


System development 


IGF1R, APP, BGN, BMP!, ESR1, LOX, ADAMFS2, FGFB1 


90304 


3.05E-02 


6 


Nucleic acid metabolic process 


APP, HNRNPUL 1, SRSF11, PFBP), ESR1, HNRNPA0 


48856 


3.18E-02 


8 


Anatomical structure development 


IGF1R, APP, BGN, BMP!, ESR1, LOX, ADAMFS2, FGFB1 


34641 


3.25E-02 


7 


Cellular nitrogen compound metabolic process 


APP, HNRNPUL1, SRSF11, PFBP), ESR1, FGFB1, HNRNPA0 


8152 


3.91 E-02 


12 


Metabolic process 


IGF1R, APP, BGN, BMP!, HNRNPUL1, SRSF11, PFBP1, ESR1, LOX, ADAMFS2, 
FGFB1, HNRNPA0 


7275 


4.02E-02 


8 


Multicellular organismal development 


IGF1R, APP, BGN, BMP!, ESR1, LOX, ADAMFS2, FGFB1 



a x represents the number of differentially expressed genes in each GO term. 
FDR, false discovery rate. 



this interaction was confirmed by immunohistochemis- 
try [22]. ESR1 can regulate bone metabolism through 
genome-wide association studies (GWAS) and it inhibits 
osteoporosis as an estrogen receptor [23]. Bioinformatics 
analysis revealed that a number of differentially ex- 
pressed genes, including ESR1 gene, are predicted to tar- 
get genes known to be important in mammalian gonadal 
development [24]. Polymorphisms at COL1A1 and 
TGFB1 and haplotypes at COL1A1 and ESR1 were found 
to be associated with bone mineral density (BMD) in a 
cohort of postmenopausal Spanish women. Moreover, 
COL1A1 polymorphisms showed significant interactions 
among them and with the VDR 3' polymorphisms [25]. 

The expression of these three genes in bone marrow 
mesenchymal cell is expected to be used in developing 
biomarkers for detecting osteoporosis and for screening 



osteoporosis risk groups. Meanwhile, the IL17RC, 
COL1A1, and ESR1 genes are related to bone protection, 
and the OPG gene is identified as balancing bone metab- 
olism by osteoclast modification. In conclusion, our re- 
sults indicate that the IL17RC, COL1A1, and ESR1 genes 
can regulate bone metabolism by the protection of the 
OPG gene in bone marrow mesenchymal cells. 

The hsa04510: focal adhesion was the most significant 
enriched pathway in the interaction network. Focal ad- 
hesions play a critical role in cell survival, migration and 
in sensing physical force. The focal adhesion pathway 
controls focal adhesion dynamics and can mediate 
reparative bone formation in vivo and osteoblast 
mechanotransduction in vitro [26]. Osteogenic differen- 
tiation is more prevalent in mesenchymal stem cells 
with a stiff, spread actin cytoskeleton and with greater 



Table 2 List of all the differentially expressed genes in the enrichment pathway 



ID 


Term 


x a 


FDR 


Genes 


hsa04510 


Focal adhesion 


21 


2.78E-09 


FLN1, ROCK2, ERBB2, IFGA11, ELK], PXN, MYL9, VEGFB, CDC42, IGF1R, RAC2, PAK2, FYN, G5K3B, ILK, 
COL6A2, COL6A1, SHC1, COL1A1, ZYX, PARVB 


hsa04722 


Neurotrophin signaling pathway 


11 


2.10E-04 


CDC42, RPS6KA2, MAP2K2, CAMK2G, GSK3B, SORF1, SHC1, MAPKAPK2, FOX03, YWHAE, ARHGDIA 


hsa04810 


Regulation of actin cytoskeleton 


12 


0.0045 


FGFR1, CDC42, ARHGEF1, PAK2, RAC2, ROCK2, MAP2K2, IFGA11, IFGB2, PXN, MYH10, MYL9 


hsa05200 


Pathways in cancer 


15 


0.007044 


FGFR1, RXRB, MAP2K2, ERBB2, SMAD3, BCL2L1, DAPK3, FGFB1, VEGFB, IGF1R, CDC42, RAC2, GSK3B, 
RARA, FRAF4 


hsa04010 


MAPK signaling pathway 


12 


0.021024 


FGFR1, CDC42, PAK2, RAC2, RP56KA2, MAP2K2, MAP3K2, JUND, ELK1, MAPKAPK2, PPP3CA, FGFB1 



a x represents the number of differentially expressed genes in each pathway. 
FDR, false discovery rate. 
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numbers of focal adhesions. Both adipogenic differenti- 
ation and chondrogenic differentiation are encouraged 
when mesenchymal stem cells have a spherical morph- 
ology associated with a dispersed actin cytoskeleton 
with few focal adhesions. Different mechanical stimuli 
can be implemented to alter these cytoskeletal patterns 
and to encourage mesenchymal stem cell differentiation 
to the desired lineage [27]. 

The bone metabolism of normal adults is in a dynamic 
equilibrium; osteoblasts synthesize new bone and osteo- 
clasts resorb old bone. If this balance is broken in vivo, 
osteoporosis is caused by insufficient bone formation 
and/or bone over-resorption [28]. Research shows that 
the main pathogenesis of osteoporosis is due to abnor- 
mal osteoblast activation and proliferation, and when 
bone absorption is more than bone formation, this nega- 
tive bone metabolism leads to osteoporosis [29]. How- 
ever, the combination of RANK and RANKL can be 
blocked by OPG, because OPG could competitive bind 
with RANKL and tumor necrosis factor-related apop- 
tosis inducing ligand (TRAIL), thus osteoclast differenti- 
ation and maturation were inhibited, and osteoclast 
apoptosis was induced. Therefore, OPG plays a key role 
against osteoporosis [30]. Osteoprotegerin, a soluble 
member of the superfamily of tumor necrosis factor re- 
ceptors, is normally secreted into marrow spaces by cells 
derived from mesenchyme. Osteoprotegerin acts as a 
decoy for osteoclast differentiation factor, which is 'both 
necessary and sufficient for osteoclast development' and 
is a critical regulator of postnatal skeletal development 
and homeostasis in humans [5]. 

The three main mechanisms of osteoporosis are an in- 
adequate peak bone mass, excessive bone resorption and 
inadequate formation of new bone during remodeling 
[31]. Therefore, to increase osteoblastic cells we can im- 
prove supplementation of osteoblastic cells, or induce 
stem cells to differentiate into osteoblastic cells for the 
treatment of osteoporosis. 

Conclusions 

By the comparison of gene chips from five osteoporosis 
patients and four normal samples of bone marrow stem 
cell, we identified genes (IL17RC, COL1A1, and £S2?i) that 
directly interact with the OPG gene. Functional and path- 
way enrichment analyses revealed that organ development 
and focal adhesion were significantly dysregulated in 
osteoporosis patients. The expressions of IL17RC, 
COL1A1, and ESR1 in bone marrow mesenchymal cell are 
expected to be used in developing biomarkers for 
detecting osteoporosis and for screening osteoporosis risk 
groups. However, further studies are still needed to con- 
firm our results because our study is based on microarray 
generated from small sample size. 
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